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Abstract 

The topological structure resulting from the network of contacts between grains 
{contact network) is studied for large samples of monosized spheres with densities 
(fraction of volume occupied by the spheres) ranging from 0.59 to 0.64. We re- 
trieve the coordinates of each bead in the pack and we calculate the average coor- 
dination number by using three different methods. We show that, in the range of 
density investigated, the coordination number is larger than 4 and it increases with 
the packing fraction. At local level we also observe a positive correlation between 
local packing fraction and number of neighbors. We discover a dependence be- 
tween the local densities of configurations with few neighbors in contact and the 
global sample-denities. This might indicate that local configurations with small 
number of neighbors are able to deform plastically when the sample is compacti- 
fying. 

1 Introduction 

The science of granular matter has a long history with classic contributions from Fara- 
day, Coulomb and Reynolds. More recently, this topic has attracted the attention of 
condensed matter physicists Q and other scientists investigating the hybrid field of 
complex systems. Indeed, the contemporary science of materials and condensed-matter 
physics has been changing in response to a new awareness of the relevance of concepts 
associated with complexity. In this context, granular materials have been sought as 
'simple' models for more complex materials and phenomena. As such, granular mate- 
rials have been typically studied either at the local level (grain scale), or as a continuous 
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Fi gure 1 1 (a) Grey-scale X-ray density map of a slice of sample C . (b) Density histogram (x-axes attenua- 
tion, y-axes frequency) for the full image volume of sample C , and (c) the same slice after phase separation 
into pore and solid phases. 

medium with specific bulk properties. The goal has always been to bridge these two 
views, in order to explain the relationships between local structure and bulk properties. 

In order to link local-structural with global, average, bulk properties it is necessary 
to develop innovative ways of analyzing the system at the grain level and to explore 
which are the possible combinations of local configurations which generate the global 
packing. Until now the empirical investigation of the geometrical structure of granular 
packings have been limited by the scarcity of accurate experimental data. Indeed, after 
the seminal works of Bernal, Mason and Scott 12] [51, it has been only very recently that 
the use of tomography has allowed scientists to 'see' the three-dimensional structure of 
such systems and explore their geometry from the grain level up to the whole packing 



In this paper we report and discuss results from an experimental investigation, by 
means of X-ray Computed Tomography, on very large samples of disorderly packed 
monosized spheres with densitie^ranging from 0.58 to 0.64. This study is the largest 
and the most accurate empirical analysis of disordered packings at the grain-scale ever 
performed |i8J. 



2 EXPERIMENTAL APPARATUS AND ANALISYS 



Detailed descriptions of the experimental methodology and apparatus are reported in 
Q m . Let us just recall that we analyzed, by means of X-ray computed tomography, 6 
samples (A, B, C, D, E, F) of monosized acrylic (polymethylmethacrylate) spheres of 
two different diameters, placed in cylindrical containers with roughened walls and with 
an inner diameter of 55 mm, filled to a height of ~ 75 mm. The fraction of volume 
occupied by the spheres in each sample and the number of spheres in the central region 
in which the analysis is performed {N) are shown in Table [T] 

'The density is defined as the fraction of the volume occupied by the balls divided by the total volume of 
the region of the space considered. 
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density 


N 


nt(1.02) 


nt(1.05) 


nt(l.lO) 


^^Deconv 


^Geom 


A 


0.586 ±0.005 


54719 


5.5 


6.7 


7.5 


5.81 


6.81 


B 


0.596 ±0.006 


15013 


5.9 


6.8 


7.7 


5.91 


7.14 


C 


0.619 ±0.005 


91984 


6.4 


7.5 


8.4 


6.77 


7.63 


D 


0.626 ±0.008 


15725 


6.0 


7.5 


8.4 


6.78 


7.73 


E 


0.630 ±0.010 


15852 


6.3 


7.6 


8.6 


6.95 


8.14 


F 


0.640 ± 0.005 


16247 


6.9 


7.9 


8.9 


6.97 


8.23 



Table 1: Sample density and their interval of variations (±) within each sample; number of 
spheres in the central region in which we perform the analysis (N); average number of neighbors 
within a given threshold distance {nt{r/d)) for different distances (r/d — 1.02, 1.05, 1.10); 
average number of neighbors resulting from the deconvolution method nneconv ; average number 
of neighbors resulting from the geometry-based method nceom- 

2.1 Reconstructed X-ray density map and material phase identifi- 
cation (segmentation) 

Figure[TJa) shows an example of an X-ray density map for a slice within a tomogram of 
sample C. The tomographic image consists of a cubic array of unsealed density values, 
each corresponding to a finite cubic portion of the sample-volume (voxel). 

The first step of the reconstruction process is to differentiate the attenuation map 
into distinct pore and solid (grain) phases. The density histogram (Fig. [ijb)) shows 
two distinct peaks associated with the two phases. The peak centered around 19500 is 
associated with the granule phase. The lower peak around 13000 is associated with the 
void phase. To distinguish between the two phases, it is sufficient to apply a threshold 
(in this case at 16000) followed by the removal of isolated solid clusters to eliminate 
noise artifacts. A comparison between the grey-scale (density map) and binarised (seg- 
mented) image of a slice of sample C is shown in Fig.[TJa) and (c). 

2.2 Determining the bead centers 

Let us consider a tomographic image with V voxels. In order to separate individual 
beads within the packing, we employ a convolution method which consists in moving 
through the (segmented) bead pack P{r) a reference sphere S{r) with radius smaller 
than the beads and measuring the overlap. The convolution of P and 5* is: 

V 

i=i 

Equation [T| can be implemented very efficiently by applying the convolution theo- 
rem which allows to transform the convolution into a product in Fourier space. 

= i{K) = T[S\ ■ T[P] (2) 

where T represents the (fast) Fourier transform. The inverse Fourier transform of Eq. [2] 
results in the beads' Image function, as expressed in Eq. [T] 
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I = T-^[i{K)] (3) 

The result of Eq.[3]is an intensity map of the overlap between the reference sphere 
and the beads, where the voxels closer to the sphere centers have higher intensities 
(in fact the peaks of the highest intensity represent the group of voxels around sphere 
centers in the original image). A threshold on the intensity map locates the groups of 
voxels surrounding the sphere centers. A visual comparison between the grains and the 
central part extracted by means of the present algorithm is shown in Fig|2] The position 
of the bead centers is then computed as the center of mass of each cluster weighted by 
the intensity of each voxel. Below, we discuss the optimal values of the threshold and 
the best choice of the reference sphere size to obtain precise estimates for the sphere 
centers' positions. 

2.3 Reference sphere size and intensity threshold 

The precision on the sphere centers' position depends on the spatial resolution of the 
tomograms and also on the size of the cluster surrounding the sphere centers after 
thresholding. The spatial resolutions of the tomograms are between one and two voxels. 
Therefore the precision on the center of mass of a cluster of v voxels must be within 
\/vX.o2/v. This suggests that, in order to minimize the error of the position of centers, 
we need to optimally choose the reference sphere size and intensity threshold so that v 
will be as large as possible. The optimal values of the parameters have been chosen by 
varying these two quantities and computing the resulting number of spheres detected 
in a given portion of the sample [8J. In Fig. [5] a reconstructed bead pack (sample E) is 
displayed. 

3 Coordination Number 

The average number of spheres in contact with any given sphere have been widely 
investigated in the literature of granular matter ||2l [5] |4] |5l [121 [Ts] [141 . Indeed, this is 




Figure 2: Beads (translusant) and the centeral part of each bead after applying the algorithm. 



4 



a very simple topological quantity which gives important information about the local 
configurations and the packing stability. 

In this paper we discuss the construction of the 'contact' networks by using three 
different methods. 

i. A thresholding method, where all the sphere centers within a given radial dis- 
tance from the centre of a sphere are considered as 'in touch'. 

ii. A deconvolution method which statistically distinguishes the contribution from 
the sphere in contact and the spheres that are close but not in touch. 

iii. A geometry-based method which employs the geometical information contained 
within the tomogram to verify if two grains are in contact or not. 

3.1 Thresholding and deconvolution methods 

From a pure geometrical point of view, the -empirical- definition of 'spheres in con- 
tact' depends on the minimum gap within which two spheres are considered in touch. 
Indeed, the exact information on the sphere positions and their shapes and sizes is em- 
pirically known only with some degree of incertitude. It is quite straightforward to 
understand that the actual number of spheres within a radial distance r from a given 
sphere (nt(r)) must increase with the radial distance: larger is the distance, and larger 
is the volume of space in which the neighboring spheres are counted. 

However, to nt{r) contribute both spheres in contact and spheres that are near but 
do not touch. It is also evident that the region concerned with the contribution from 
the spheres in contact can be only within a small region with size compatible with the 
sphere polydisperity. The 6 samples analyzed in |i8J have polydispersities around 2% 
of their diameters, therefore it is reasonable to assume that above the radial distance 




Figure 3: The network of connected beads in a reconstruction of sample E. 
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1.02d all the contributions to nt{r) come only from spheres not in contact. We can 
therefore (over-)estimate the number of grain in 'contact' by counting all couple of 
spheres with centers within radial distances r < 1.02d. The results for the average 
numbers of spheres in 'contact' at threshold distances r/d — 1.02, 1.05 and 1.10 are 
reported in Table [T] 

The thresholding procedure has the defect to count both the contributions from the 
spheres in contact and from other spheres which are near but do not touch. The chal- 
lenge is to distinguish between these two different contributions. In |8 1 a deconvolution 
method was introduced to overcome this problem. THe results of such an analysis are 
reported in the column noeconv in Table [T] 

3.2 Geometry-based method 

Differently from the methods discussed in the previous section, this method detects 
the touching spheres by using directly the information from the segmented image of 
granular packs. Indeed, the segmented image provides a full 3D structural information 
of the bead samples and therefore one can directly determine whether two beads are 
in contact or not. However, this method, depends highly on the resolution of the the 
tomogram. Since there is an uncertaint}|^of 1 to 2 voxels in determining the grain-grain 
boundary, we expect to detect neighbors within radial distances of 1 — 2 voxels^/^^. This 
means that for beads with diameters of 25 voxels^/'^ (B, D, E, F) |[8|, the uncertainty is 
within radial distances of 1.04 — 1.08 x d. Whereas an uncertainties of 1.03 — 1.07 x d 
are expected in samples (A, C) where the diameter is 30.81 voxels^ |8 1. 

The results for the average number of neighbors in contact, computed by means 
of this geometrical-based method are reported in Table. [T] (nGeom)- An example of 
contact network resulting from this analysis is shown in Fig. [3] A comparison between 
these data and the results from the deconvolution method shows a very good agreement 
between these two techniques when the implicit thresholding due to the voxel-precision 
is considered. As expected, we find that the geometry-based method overestimates the 
number of actual contacts giving results equivalent to the tresholding methods when 
neighbors within radial distances between 1.05 and 1.07 bead diameters are counted. 

In Fig |4] the values of the average number of contacts are reported as a function of 
the sample densities. From the figure it is clear that all the methods reported in this 
paper and also other independent experimental investigations, give average numbers of 
contacts larger than 4 which increase with the sample-density. Such dependence on 
the packing density (Fig|4|i has important conceptual impUcations especially for what 
concerns the kind of mechanical equihbrium in the system fS",^. 

4 Local density fluctuations and local number of neigh- 
bors 

In the previous section we have observed that the average number of neighbors in 
contact is affected by the sample-density. Here we look at the same relation but at local 

^This is due to the diffraction of X-ray at tlie boundary 
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Figure 4: Number of neighbors in contact vs sample's density. The dotted line is an 
extrapolation to ric = 4 at p = 0.55. 



level. To this end let us first define the local density as the fraction of Voronoi cell 
occupied by each sphere in the sample. Such a quantity varies from sphere to sphere. 
The distribution of local densities are reported in Fig|5]^right). One can observe that, 
there are no local configurations denser than 0.75 and there are very few configurations 
with density lower than 0.5. The distributions have bell-like shapes with averages 
around the sample-densities and standard deviations between 5 and 6 % (decreasing 
with increasing sample-density). The number of neighbors also varies from sphere to 
sphere. Figure |5jleft) reports such a distribution for the 6 samples A-F. As one can see, 
the distributions of the number of neighbors are also bell-shaped around the sample- 
averages and a dependence on the sample-density is observed with denser samples 
showing larger connectivities. 

Here we are interested in investigating the relation between the number of neigh- 
bors per sphere and the local densities. Such a relation is reported in the two figures 
Fig. [6|a),(b) where both the average number of neighbors v.s. local density and the 
average local density v.s. number of neighbors are reported. Let us first point out that, 
these two plots are not simply the reciprocal inverse of each other Indeed they corre- 
spond to two opposite gatherings. The Fig. |6|a) is obtained by counting the number 
of neighbors (within a threshold distance of 1.02(i) in local configurations with a given 
density. Conversely, Fig. [6|b) reports the average local densities of configurations with 
given number of neighbors (1.02d threshold distance criteria). In both cases the error 
bars correspond to the standard deviations of the respective distributions. 
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Figure 5: (left) Distribution of the number of neighbors per sphere (threshold method 
with distance 1.02d). (right) Distribution of local density (fraction of occupied Voronoi' 
cells). 




Figure 6: (a) Average number of spheres counted around local configurations with 
given density, (b) Average local density of configurations with a given number of 
neighbors. In both figures the error-bars are the standard deviations. 



From these figures it is evident that there is a positive correlation between the num- 
ber of neighbors and the density of local configurations. This is not surprising, close 
packed configurations must have a larger number of neighbors within a given distance 
than loose packed configurations. Indeed, it is known that the most compact configu- 
ration can have density 0.754 with 12 neighbors in contact. Whereas, on the opposite 
limit. Fig. |6|b) shows that there can be some configurations with only two neighbors 
and local densities in a range between 0.5 to 0.62. Beyond such rather intuitive out- 
comes, two interesting fact emerge from Figs. |6] 

1) local configurations with larger density have a larger number of neighbors, but 
this is essentially independent on the global sample density (see Fig. |6]^a)); 
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2) configurations with larger number of neighbors also show larger local densities 
but in this case a global dependence on the sample density is also observed (see 

Fig.igb)). 

Such a relation between local and global quantities is quite remarkable. These data 
show that local configurations with a small number of neighbors have larger, or smaller, 
local densities depending on the global density of the sample. This might suggest that 
the geometrical structure of sphere packings is adapting 'plastically' by expanding or 
contracting some local configurations while constraining mechanical stability. 

5 Conclusion 

The study of the topological structure of the contact network at the grain level shows 
that the average number of spheres in contact increases with packing density and lies 
between 5.5 and 7.5 in the range of densities between 0.58 and 0.64. An extrapolation 
to the random loose packing density (p — 0.55) suggests that at this density the sys- 
tem could have an average number of 4 neighbors per sphere. This might imply the 
possibility of a rigidity percolation transition taking place at the random loose packing 
limit lH). It should be however considered that the current methods to estimate the 
number of contacts have an excessively large incertitude that does not allow to validate 
-or invalidate- the extrapolation to = 4 at p = 0.55. Further studies are needed to 
clarify this point. One way to reduce such incertitude on the geometrical number of 
neighbors is to consider in touch only grains that exchange some stress at the interface. 
This can be done both by means of numerical simulations or by direct measurements 
of stress in three-dimensional deformable grains. But, this will be the object of future 
research. 

At local level the we observed that there is a positive correlation between the num- 
ber of neighbors and density: spheres with a larger number of neighbors occupy a 
larger fraction of their local Voronoi cells. We also observed that configurations with 
small number of neighbors (2-6) are more sensitive to the global packing density. This 
suggests that, while constraining global mechanical stability, the system is able to com- 
pactify by squeezing 'plastically' the configurations with smaller number of neighbors. 
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